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Abstract.  Time  accurate  numerical  simulations  of  a  high  temperature  source  cloud  of  gas  expanding  into  an  ambient 
atmosphere  are  performed  using  a  multiple  temperature  gas  model  and  the  direct  simulation  Monte  Carlo  (DSMC) 
method.  The  multi-temperature  approach  uses  continuum  conservation  equations  derived  from  the  Boltzmann  equation 
via  first-order  Chapman-Enskog  expansion  and  zero-order  non-isotropic  velocity  distribution  function.  These  equations 
are  solved  numerically  using  a  kinetic  flux  splitting  method  for  inviscid  fluxes  and  a  central  difference  scheme  for  the 
viscous  fluxes  in  a  time  accurate  manner.  The  DSMC  technique  is  a  well-established  particle  approach  for  rarefied  flows. 
The  source  and  ambient  densities  and  temperatures  are  varied  to  establish  the  range  of  applicability  and  relative  accuracy 
of  the  two  methods.  For  most  of  the  source  and  the  ambient  parameter  considered,  the  two  methods  give  similar 
predictions  of  the  basic  flowfield  evolution,  and  of  the  degree  of  translational  non-equilibrium  that  arises. 


INTRODUCTION 

The  expansion  of  a  highly  energetic  gas  cloud  into  a  rarefied  environment  is  a  phenomenon  that  arises 
in  many  practical  problems,  including  supernova  explosions,  high  altitude  rocket  plumes,  and  active  space 
experiments  where  interaction  of  the  artificial  plasma  and  neutral  gas  cloud  with  the  ambient  atmosphere  is 
studied.  During  the  rapid  expansion  of  the  gas  cloud,  the  spherical  gas  front  or  contact  surface  propagates 
outward  through  the  ambient  atmosphere,  while  a  rarefaction  wave  propagates  inward  into  the  source  cloud. 
Recent  time-accurate  continuum  simulations1  of  a  high-temperature  metallic  gas  cloud  at  high  altitude  with 
radiation  transport  show  the  flow  field  generally  includes  continuum,  transitional,  and  rarefied  regions 
depending  on  the  absolute  and  relative  densities  and  thermal  energies  of  the  source  cloud  and  ambient 
atmosphere.  Thus  it  seems  that  particle  simulation  methods  such  as  the  direct  simulation  Monte  Carlo 
(DSMC)  method2  would  be  more  suitable  for  simulating  these  flows.  However,  the  unsteady  nature  of  these 
flows  makes  the  particle  approach  a  formidable  computational  task.  Other  features  that  add  further 
complexity  to  the  particle  treatment  are  the  presence  of  very  dense  gas  in  the  source  region  and  the  potential 
for  phase  change  as  well  as  coupled  radiative  transport.  Moreover,  the  dilute  gas  assumptions  made  in  these 
particle  simulation  methods  make  them  unsuitable  for  dense  regions  of  the  cloud.  This  suggests  that  at 
present  a  continuum  method,  which  to  some  extent  accounts  for  the  non-equilibrium  nature  of  the  flow,  may 
provide  efficient  and  accurate  simulations  of  the  gas  cloud  expansion. 

The  purpose  of  this  study  is  to  improve  the  fundamental  understanding  of  the  initial  evolution  of  the 
gas  cloud  expansion,  including  the  formation  of  the  shock  wave  in  the  ambient  atmosphere  and  the  degree  of 
translational  non-equilibrium  that  arises.  Time  accurate  numerical  simulations  of  a  rapidly  expanding  high 
temperature  gas  cloud  are  performed  using  both  a  multi-temperature  gas  dynamics  model3"4  and  the  DSMC 
method2.  The  influence  of  the  source  radius  and  source  and  ambient  density  are  evaluated.  Comparison  of 
results  for  the  two  methods  establishes  the  range  of  applicability  and  the  relative  accuracy  of  each. 
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COMPUTATIONAL  APPROACHES 


Conservation  equations  of  the  multi-temperature  gas  dynamic  model  in  spherical  symmetric  coordinates 
(Ref.4)  are  solved  numerically  in  a  time  accurate  manner.  The  inviscid  part  of  the  flux  is  calculated  at  each 
cell  interface  using  a  kinetic  flux  splitting  method. 1  A  local  ellipsoidal  distribution  function  of  the  form 
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is  then  defined  for  each  cell  and  used  to  formulate  analytical  expressions  for  the  inviscid  fluxes  in  the  radial 
(r)  and  perpendicular  (p)  directions.  The  c  c and  c  are  the  components  of  the  random  molecular  velocity 
in  the  above  expressions.  As  with  other  flux-splitting  techniques,  this  inviscid  flux  is  formally  first-order 
accurate  in  space.  The  second-order  accuracy  in  space  is  achieved  by  using  a  MUSCL  variable-extrapolation 
scheme5  in  conjunction  with  a  minmod  slope  limiter  to  ensure  monotonicity.  The  viscous  fluxes  are 
calculated  with  a  second-order  accurate  central  difference  scheme.  Fourth-order  accuracy  in  time  is  obtained 
using  a  Runge-Kutta  time  integration  method6.  The  time  step  chosen  is  the  minimum  of  a  Courant-based 
time  step  and  a  step  based  on  the  magnitude  of  energy  exchange  between  the  radial  and  perpendicular 
directions. 

The  DSMC  calculations  are  performed  using  a  conventional  one-dimensional  method  well  described  in 
reference  2.  Some  features  of  the  current  implementation  are  given  in,  e.g.,  reference  7.  Profiles  of 
macroscopic  properties  at  a  given  time  are  generated  as  instantaneous  snapshots  of  sampled  molecular 
quantities. 


RESULTS  AND  DISCUSSION 

The  source  cloud  is  assumed  to  be  spherical  with  uniform  initial  density  and  temperature.  The  source  is 
pure  argon  with  an  assumed  Maxwellian  (inverse  fourth  power  law)  inter-atomic  potential 

(  Tref  =  273 K and  pret  =  2.5 1 7 x  10  5  kg! ms).  Two  values  of  source  gas  density  ps  at  a  temperature  of  Ts  = 

3000  K  are  considered  in  the  present  calculations.  The  outer  boundary  condition  is  assumed  to  be  the 
standard  atmosphere  at  the  given  altitude  (70  or  85  km).  To  focus  on  cloud  expansion  and  not  on  molecular 
diffusion  and  other  multi-species  issues,  the  ambient  atmosphere  is  simulated  as  pure  argon  at  the  ambient 
atmospheric  density  and  temperature. 


In  the  multi-temperature  continuum  calculations  the  initial  spherical  source  volume  and  flow  domain 
are  divided  into  uniformly  distributed  spherical  cells  of  constant  width  0.0002m.  One  cell  is  added  to  the 
outer  boundary  after  each  time  step,  allowing  the  grid  to  grow  in  time  as  the  gas  cloud  expands.  This 
procedure  reduces  the  computational  task  during  early  phases  of  the  solution  because  computation  may  need 
very  large  number  of  cells  at  later  times. 
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The  DSMC  calculations  use  a  uniform  radial  grid  with  cell  dimension  comparable  to  the  initial  source 
mean  free  path.  Collision  sub  cells  were  utilized  to  further  improve  the  spatial  fidelity.  Maxwell  molecule 
collision  interaction  was  approximated  using  the  Maxwell  limit  of  the  VHS  model.  Sensitivity  studies  (not 
shown)  confirm  the  results  to  be  grid  and  timestep  independent.  Since  spatial  weight  factors  are  not  used, 
initial  molecular  population  of  a  cell  is  controlled  only  by  the  local  gas  number  density  and  the  cell  volume. 
Control  of  statistical  noise  present  in  the  sampled  snapshots  is  possible  in  this  case  only  by  recourse  to  large 
numbers  of  simulated  molecules.  The  typical  simulation  required  40  million  molecules,  1000  cells,  10 
subcells  per  cell,  and  2500  timesteps. 


Figures  1-4  show  the  evolution  of  the  density,  temperature,  and  velocity  profiles  for  a  source  cloud 
expanding  into  a  70  km  altitude  ambient..  The  initial  source  density  and  radius  for  this  case  are  ps  =  0.005 
kg/m3  and  r,  =0.1  m,  corresponding  to  a  source  Knudsen  number  of  Kns  ~  5.6xl0'4.  The  source  to  ambient 
density  ratio  is  approximately  60.  The  initial  source  energy  is  19.6  J  with  energy-density  of  4.68xl03  J/m3. 
These  figures  show  that  a  well-defined  shock  wave  propagates  outward  in  the  ambient  atmosphere.  The 
shock  is  followed  by  a  contact  surface,  which  separates  the  expanding  gas  from  the  ambient  environment.  A 
strong  rarefaction  wave  travels  behind  the  contact  surface  towards  the  center  of  the  source.  The  final  time 
shown  (1=1x1  O'4  s)  corresponds  to  approximately  one  acoustic  time  (based  on  the  initial  source  radius  and 
speed  of  sound). 

Translational  non-equilibrium  arises  in  the  shock  wave  as  the  gas  cloud  expands  (Fig.  3).  The  ratio  of 
the  radial  temperature  to  the  parallel  temperature  (Tr  /Ts)  as  a  function  of  normalized  gas  cloud  radius  is 
shown  in  Fig.  3  for  three  different  times.  The  maximum  value  of  this  ratio  is  more  than  two  which  occurs  in 
the  shock  wave  region  (Fig.  3).  This  suggests  that  the  shock  wave  region  is  in  a  highly  non-equilibrium 
thermal  state.  The  majority  of  the  expanding  gas  cloud  remains  in  translational  equilibrium,  due  to  the 
relatively  high  source  density.  A  slight  overshoot  is  observed  in  the  radial  temperature  just  behind  the  shock 
wave  (Fig.  2),  but  no  overshoot  is  present  in  the  perpendicular  temperature  (not  shown). 

The  radial  velocity  increases  linearly  with  radius  through  the  cloud  for  most  of  the  flow  field  (Fig.  4). 
Also,  the  radial  velocity  gradient  decreases  with  time  as  the  gas  cloud  expands  into  the  ambient  atmosphere. 


Altitude  70  km 


Altitude  70  km 


Fig.  1  Density  Profiles 
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Altitude  70  km 


Fig.  3  Parallel  Temperature  Profiles 


Fig.  4  Radial  Velocity  Profiles 


DSMC  calculations  predict  a  slightly  thicker  shock  than  predicted  by  the  multi-temperature  continuum 
model  (Figs.  1-4).  The  DSMC  results  are  expected  to  be  more  accurate  than  the  continuum  method  in  the 
highly  non-equilibrium  leading  edge  of  the  shock.  Also,  the  difference  in  shock  wave  thickness  is  slightly 
more  pronounced  in  the  temperature  profiles  because  temperature  is  calculated  based  on  higher  moments  of 
the  distribution  function.  The  DSMC  results  show  relatively  large  scatter  near  the  origin  due  to  the  small 
cell  volumes  present,  but  this  does  not  obscure  interpretation  of  the  remainder  of  the  flowfield.  The  overall 
qualitative  and  quantitative  comparison  of  the  flow  field  predictions  by  the  two  methods  is  very  good  for 
this  case. 

A  second  set  of  calculations  (Figs.  5-8)  is  also  performed  at  70  km  altitude  for  expansion  of  a  gas  cloud 
with  five  times  smaller  source  density  and  half  the  radius  of  the  previous  case.  For  this  case,  rs  =  0.05  m  and 
ps  =  0.001  kg/m3,  resulting  in  Kns  ~  5.6xl0"3  and  p,  lpa  ~  12.  The  source  energy  is  0.489  J,  and  the  energy 
density  is  934  J/m3.  These  calculations  illustrate  the  effect  of  initial  source  energy  on  the  gas  cloud 
expansion  for  given  ambient  conditions.  Figures  5-8  show  the  same  general  features  of  the  gas  cloud 
expansion  as  are  observed  in  the  first  set  of  calculations  (Fig.  1-4)  except  that  a  shock  wave  does  not  form 
in  the  ambient  because  there  is  not  enough  source  energy  to  raise  the  expansion  front  velocity  to  high 
supersonic  value  (Fig.  8).  Comparison  of  Fig.  6  and  Fig.  7  of  the  temperature  profiles  shows  considerable 
thermodynamic  non-equilibrium  in  the  expansion  front  region. 


Altitude  70  km 


Altitude  70  km 


Fig.  6  Radial  Temperature  Profiles 
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Altitude  70  km 


Altitude  70  km 


Fig.  7  Parallel  Temperature  Profiles 


In  order  to  see  the  effect  of  ambient  rarefaction  on  the  gas  cloud  expansion,  a  third  set  of  calculations  is 
performed  for  the  source  conditions  used  in  the  first  set  (Fig.  1-4),  but  with  an  ambient  corresponding  to  an 
altitude  of  approximately  85  km.  The  ambient  values  of  density  and  temperature  are  pa  =  8.2196  x  10"6  kg/m3 
and  Ta  =188.893  K,  giving  a  density  ratio  p,  tpa  ~  600.  Figures  9  and  10  show  normalized  plots  of  density 
and  radial  temperature  as  a  function  of  normalized  radial  distance.  Fig.  9  shows  that  the  multi-temperature 
model  predicts  a  well-defined  thin  shock  in  the  ambient,  whereas  the  DSMC  calculations  predict  a  thick  and 
diffused  shock  wave.  The  multi-temperature  model  cannot  resolve  the  large  non-equilibrium  effects  present 
in  the  leading  edge  of  the  disturbance,  and  for  this  reason  the  radial  temperature  in  the  shock  region 
predicted  by  the  multi-temp  model  is  lower  than  predicted  by  the  DSMC  method  (Fig.  10).  The  difference  in 
parallel  temperatures  in  the  shock  wave  region  (not  shown)  is  much  smaller  than  for  radial  temperatures. 
Also,  the  difference  in  radial  velocity  predicted  by  the  two  methods  in  the  shock  wave  region  (not  shown)  is 
small.  The  overall  comparison  of  flow  field  predictions  by  two  methods  is  reasonably  good. 


Altitude  85  km  Altitude  85  km 


Fig.  9  Density  Profiles 


Fig.  10  Radial  Temperature  Profiles 
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CONCLUSIONS 


Computations  carried  out  for  a  variety  of  source  and  ambient  conditions  demonstrate  the  effects  of 
source  energy  and  ambient  rarefaction  on  the  evolution  of  the  gas  cloud.  For  the  low  source  energy 
conditions  at  70  km  altitude,  no  shock  wave  propagates  in  the  ambient  because  the  gas  front  velocity  is 
barely  supersonic.  Higher  source  energy  conditions  at  the  same  altitude  result  in  a  well  defined  shock  wave, 
with  considerable  translational  non-equilibrium  present  in  the  shock  wave  region.  The  multi-temperature 
continuum  model  and  the  DSMC  method  predict  a  radial  temperature  overshoot  in  the  gas  front  region  for 
the  dense  and  high  energy  source  conditions  at  70  km  altitude.  No  radial  temperature  overshoot  is  predicted 
for  the  same  source  expanding  at  85km  altitude  where  ambient  atmosphere  is  more  rarefied  by  an  order  of 
magnitude.  The  comparison  of  flow  structure  of  the  gas  cloud  expansion  predicted  by  the  multi-temperature 
model  and  the  DSMC  method  shows  the  multi-temperature  model  is  a  reasonably  good  computational  tool 
for  low  to  moderate  rarefied  and  non-equilibrium  conditions. 


REFERENCES 

1.  Dogra,  V.  K.,  Nance,  R.  P.,  Taylor,  J.  C.,  Swaminathan,  P.  K.,  Erlandson,  R.  E.,  and  Meng,  C.-I.,  “Modeling  of 
Gas  Cloud  Expansion  at  High  altitude  with  Radiation  Transport,”  Journal  of  Thermophysics  and  Heat  Transfer, 
Vol.  14,  No.  3,  July-September  2000. 

2.  Bird,  G.A.,  Molecular  Gas  Dynamics  and  the  Direct  Simulation  of  Gas  Flows,  Clarendon  Press,  Oxford, 
England,  1994. 

3.  Candler,  G.  V.,  Nijhawan,  S.,  Bose  Deepak,  and  Boyd,  I.D.,  “A  Multiple  Translational  Temperature  Gas 
Dynamics  Model,”  Physics  of  Fluids,  Vol.  6,  No.  11,  1994. 

4.  Dogra,  V.  K.,  et.,  ah,  “Simulations  of  Gas  Cloud  Expansion  Using  a  Multi-Temperature  Gas  Dynamics  Model,” 
RGD  22nd  International  Symposium,  AIP  conference  proceedings,  vol.  585,  (2000). 

5.  Van  Leer,  B.,  “  Towards  the  Ultimate  Conservative  Difference  Scheme.  A  Second  Order  Sequel  to  Godunov’s 
method,”  Journal  of  Computational  Physics,  Vol.  32,  1979,  pp.  263-275 

6.  Jameson,  A.,  Schmidt,  W.,  and  Turkel,  E.,  “Numerical  Simulations  of  the  Euler  Equations  by  Finite  Volume 
Methods  using  Runge-Kutta  Time  Stepping  Schemes,”  AIAA  Paper  81-1259,  1981. 

7.  Wysong,  I.  J.  and  Wadsworth,  D.  C.,  “Assessment  of  direct  simulation  Monte  Carlo  phenomenological 
rotational  relaxation  models”,  Physics  of  Fluids,  Vol.  10,  Number  11,  1998,  pp.  2983-2994. 


412 


